Convergence condition (not converge for every matrix); Slow for exact solutions
Tensor Calculus
Basic Concepts
1st-Order Derivatives
If , then
If , then
2nd-Order Derivatives
If , then (Hessian is symmetric, tangent stiffness)
Taylor Expansion
If , then
If (vector func) , then
For is s.p.d., second order derivative > 0 => interesting properties (to be discussed)
Lecture 3 Rigid Body Dynamics
(Single rigid body: dynamics / rotation / …)
Rigid bodies: Assume no deformations
The goal of simulation is to update the state var. ove
Translation Motion
For translation motion, the state variable contains the position and the velocity ( - Mass; Force can be the function of pos, vel, t, …) -> Solve integral
Integration Methods Explained
By def, the integral of is the area.
Explicit Euler (1st-order accurate) sets the height at
Use Taylor Expansion: (Error: )
Implicit Euler (1st-order accurate): sets the height at
Taylor: (Error: )
Mid-point (2nd-order accurate): sets at
Taylor:
Final Method in this case: Semi-implicit (Mid-point)
Velocity -> Explicit; Position -> Implicit
Alternative: Leapfrog Integration
v and x not overlap (Mid-points)
Type of Forces
Gravity Force:
Drag Force: ( - drag coefficient) -> Reduced by following
Use a coefficient to replace the drag force: ( - decay coefficient)
Translation Only Simulation
Steps: (Mass and Timestep are user spec var)
Rotation Motion
Rotation Representations
Rotation Represented by Matrix
Suitable in graphics, rotation for vertices
Not suitable for dynamics:
Too much redundancy: 9 elem, 3 dof
Not intuitive
Defining its time derivative (rotational vel) is difficult
Rotation Represented by Euler Angles
Use 3 axial rotations to represent one general rotation. Each axial rotation uses an angle.
Used in Unity. (the order is rotaion-by-Z / X / Y) Intuitive.
Not suitable for dynamics:
Lose DOFs in certain statues: Gimbal lock
Defining its time derivative is difficult
Rotation Represented by Quaternion
Complex multiplications: In the complex system, two numbers represent a 2D point. => Quaternion: i, j, k are imaginary numbers (3D space) Four numbers represent a 3D point (with multiplication and division).
Arithematic
Let be a quaternion made of 2 parts: a scalar and a 3D vector for
Scalar-quaternion Multiplication
Addition/Subtraction (Same as vector)
Multiplication
Magnitude
In Unity: provide multiplication, but no addition/subtraction/…; Use w, x, y, z -> s, v
Representation
Rotate around by angle
Convertible to the matrix:
Rotation Representations in Unity
Rotation Motion
Use a 3D vector to denote angular velocity:
The dir of -> the axis; The magnitude of -> the speed (sim to the representation of quaternion)
Torque and Inertia
Torque
(Original state: -> Rotated: , is a force)
The rotational equiv of a force, describing the rotational tendency caused by a force.
: perpendicular to both and ; proportional to and ; proportional to (the angle of the two vectors)
Inertia
Describes the resistance to rotational tendency caused by torque (not const)
Left side (heavier point far away from the torque) has higher resistance (inertia) to the rotational tendency, slower rotation
Ref state inertia, change with rotation (dep on pose). But no need to re-compute every time
( - 3x3 identity matrix)
Use torque to represent
Torque on a spec point:
Total torque:
The rotational equivalent of mass is called inertia (Result: 3x3):
Reference inertia:
Current inertia:
Rigid Body Simulation
Translational and Rotational Motion
Translation (Linear)
States: velocity and position (transform.position in Unity)
Physical Quantities: mass and force
Rotational (Angular): Better normalize when (in Unity automatically)
States: angular velocity and quaternion (transform.rotation in Unity)
Physical Quantites: inertia and torque
Rigid Body Simulation Process
In Unity: No 3x3 matrices, only 4x4 (use 4x4 and set the last col / line ); Provide .inverse to inverse; …
Implementation
In practice, we update the same state var
Issues
Translation is easier, code translation first
Using a const first while testing update , in this case the object will spin constantly
Gravity does NOT cause torque (except for air drag force)
Lecture 4 Rigid Body Contacts
Particle Collision Detection and Response
Distance Functions
Signed Distance Function
Use a signed distance func to define the distance indicating which side as well (corresponding to 0 surface)
Examples
Intersection of Signed Distance Functions
=> Bool operations
if and and then inside// all val are negative, the max one is the closestelse// not relavent (no collision)
Union of Signed Distance Functions
if or then inside// approximate -> correct near outer boundaryelse outside
-> We can consider collision detection with the union of two objects as collision detection with two separate objects
Penalty methods
(Implicit integration is better)
Quadratic Penalty Method
Check if collide - Yes -> Apply a force at the point (in the next update)
For quadratic penalty (strength) potential, the force is linear
Problem: Already inside -> cause artifacts
=> Add buffer help less the penetration issue (cannot strictly prevent)
Problem: too low -> buffer not work well; too high -> too much force generated (overshooting)
Log-Barrier Penalty Method
Ensures that the force can be large enough, but assumes never happens => By adjusting
Always apply the penalty force: ( - Barrier strength)
Problems: cannot prevent overshooting when very close to the surface; when penatration occurs, the penatration will be higher and higher (-> smaller step size -> higher costs)
=> Log-Barrier limited within a buffer as well to solve
Frictional contacts are difficult to handle
Impulse method
Update the vel and pos as the collision occurs
Changing the position:
Changing the velocity: ( - bounce coefficient, , - frictional decay of vel)
should be minimized but not violating Coulomb’s law
Therefore:
Can precisely control the friction effects
Rigid Collision Detection and Response by Impulse
Rigid Body Collision Detection
When the body is made of many vertices, test each vertex: (from the mass center to the vertices)
=> detection: transverse every point if
Rigid Body Collision Response by Impulse
Vertex : ( - linear vel; - angular vel)
But cannot modify and directly since they are not state var. (in Lec.3: 4 var are vel of mass center / pos of mass center / rotational pose / angular vel)
Applying an impulse at vertex : ( (Newton’s Law); - torque induced by )
Cross Product as a Matrix Product
Convert the cross prod into a matrix prod
In our case:
Therefore, replace with some matrix . Finally can be computed with the following equations.
Implementation
Other details:
For several vertices in collision, use their average
Can decrease the restitution to reduce oscillation
Don’t update the position: not linear
Shape Matching
Basic Idea
Allow each vertex to have its own velocity, so it can move by itself
Move vertices independently by its velocity, with collision and friction being handled (use the impulse method for every point)
Enforce the rigidity constraint to become a rigid body again (IMPORTANT)
Mathematical Formulation
Want to find and ( - mass center): want the final rigid body (a square) close enough to the trapozoid
( - a matrix, not only corresponding to rotation, since the mass center was set to 0; - The objective, )
For mass center and matrix (Find derivatives):
Remind: Singular value decomposition: (rotaion, scaling and rotation)
Rotate the object back before the final rotation: (Local deformation: )
Implementation
Physical quantities are attached to each vertex, not to the entire body.
(The function of is provided, ultilizing the polar deposition technique)
Properties:
Easy to implement and compatible with other nodal systems: cloth, soft bodies, particle fluids, …
Difficult to strictly enforce friction and other goals. The rigidification process will destroy them (may require iter)
More suitable when the friction accuracy is unimportant, i.e., buttons on clothes
Lecture 5 Physics-Based Cloth Simulation
A Mass-Spring System
Spring Systems
An Ideal Spring
Satisfies Hooke’s Law: The spring force tries to restore the rest length ( - Spring Stiffness)
Multiple Springs
The energies and forces can be simply summed up
Structures in Simulations
Structured Spring Networks
Unstructured Spring Networks
Unstructured triangle mesh
-> the edges into spring networks (usually in cloth simulations)
Blue lines for bending resistance (every neighboring triangle pair)
(Notations: - Mass of vertex ; - Edge list; - Edge length list (pre-computed))
For every vertex:
To compute Spring Forces: (For every spring )
(Spring index and vertex index )
Problem
Overshooting: when and/or is too large
Naive solution: Reduce => Slow down the simulation
Implicit Integration
General Scheme (Euler Method)
Integrate both and implicitly ( - Mass matrix (usually diagonal))
& are unknown for the current time step => Find the x and v:
Assume dep only on (homonomic): Solve the eqn (Problem: may not be a linear function)
The equation equiv to the following (where ) => Optimization prob => Numerical schemes (Usually only conservative forces can use this energy)
Because: (applicable for every system; The first order der of reaches 0 for the min pt.)
The Optimization Problem
Newton-Raphson Method
Solving optimization problem: ( is Lipschitz continuous)
Given a current we approx the goal by (Taylor Expansion)
(For 2D: )
Steps:
Initialize
For
(For 2D: )
If is small then break
Newton’s Method finds extremum, but it can be min or max => finds 2nd order derivative (at local min: ; at max: )
For a function which second order derivative always larger than 0 ( everywhere) => has only one minimum
Simulation by Newton’s Method
For simulation:
Derivative:
Force Hessian Matrix: ( - Energy Hessian)
Steps:
Initialize , often as or
For :
Solve
If is smallthen break
Spring Hessian
Hessian matrix is a second order derivative (sim to ) => if p.d. so that the ONLY min point is found and has no maximum (sufficient but not neccesary condition)
The matrix in the last part is 3N x 3N, every vertex is a 3D vector. The first is at , the in the first line is at , …
Positive Definite: Dep on every (3 x 3): @ Lec.2, P48
where & are already s.p.d. (Proof by multiplying by and at both sides)
But can be negative when (when compressed)
Conclusion: when stretched is s.p.d. and when compressed may not be s.p.d. => may not be s.p.d. either
(for smaller => more p.d., actually sim to explicit integration with smaller time step => more stable)
Positive Definiteness of Hessian
When a spring is compressed, the spring Hessian may not be positive definite => Multiple minima
Only in 2D/3D: In 1D: and
Enforcement of P.D.
Some linear solver may require must be p.d. in
Solution:
Drop the ending term when :
Choi and Ko. 2002. Stable But Responive Cloth. TOG (SIGGRAPH)
Linear Solvers
Solving
Jocobi Method
For :
// Residual error
if then break // Convergence condition
// Update by , the diagonal of
vanilla Jocobi method () has a tight convergence req on : Diagonal Dominant
Baraff and Witkin. 1998. Large Step in Cloth Simulation. SIGGRAPH.
One of the first papers using implicit integration. The paper proposes to use only one Newton iteration, i.e., solving only one linear system. This practice is fast, but can fail to converge.
Bending and Locking Issues
The Bending Spring Issue
A bending spring offers little resistance when cloth is nearly planar, since its length barely changes.
A Dihedral Angle Model
Every vertex of the 4 will be under a force as a function of :
and should be in the normal directions &
Bending doesn’t stretch the edge should be orthogonal to the edge (in the span of & )
Newton: , means & are in the span of &
Conclusion:
; (Cross prod tells the dir. (for normal should be normalized) )
Force:
Planar Case: (The magnitude of the cross product tells the area ())
Non-planar Case: (The initial angle is at rest, planar case )
Explicit Derivative is difficult to complete / No energy
A Quadratic Bending Model
2 assumptions: planar case (OK for cloth); little stretching (mainly bending)
Energy function: (vector * matrix * vectorT) ( is 3x3 identity)
Actually finding the laplacian (/ curl), when flat, no curl =>
is a constant matrix => The function is quadratic
Pros & Cons
Easy to implement (even in implicit method)
No longer valid if stretches much; Not suitable if the rest config is not planar -> (projective dynamics model / cubic shell model)
The Locking Issue
In the mass-spring / other bending models, assuming cloth planar deformation and cloth bending deformation are independent
In zero bending case: LHS - OK; RHS - For stiff spring NO (Locking Issue) => short of DoFs.
For a manifold mesh, Euler formula: #edges = 3#vertices - 3 - #boundary_edges
If edges are all hard constraints: DoFs = 3 + #boundary_edges
=> no perfect solution
Shape Matching
..
Lecture 6 Constraint Approaches: PBD / PD / …
Strain Limiting and Position Based Dynamics
The Stiffness Issue
Real-world fabric resist strongly to stretching
But in simulation, only increasing the stiffness can cause: (More expensive / time-costing)
Explicit integrators - Unstable
-> Smaller timesteps and more computational time
Linear systems invloved in implicit integrators will be ill-conditioned
-> More iter and computational time
A Single Spring
If a spring is infinitely stiff -> Length = const
Def a constraint func: (Rest length )
Def a proj func: suppose in a space for a pos of want to move into a rational area (in blue, with boundary ) with a shortest path: want the projection
The 2 new points’ substraction equals the original length => satisfying the constraint
By default , but can also set for stationary nodes (can be just ignored)
Multiple Springs
Gauss-Seidel Approach
Approach every spring sequentially in a certain order (needs a lot of iter to converge)
Cannot ensure the satisafction of every constraint. More iter give much closer results (to the constraint)
More similar to the stochastic gradient descent (in ML, with spec order)
The order mattrers => cause bias and affect convergence
Jocabi Approach
Projects all edges simultaneously (good for parallization) and then linearly blend the results
sum up all and update together with a weighted average.
Lower convergence rate
More iter give better results
Position Based Dynamics (PBD)
Based on the proj func, similar to shape matching
& update as simple particle system (similar to Shape matching in rigid body collision)
PBD: Add constraints
: use Gauss-Seidel / Jocobi /… to find the projection
: use new to find
Properties:
The stiffness behavior is subject to non-physical factors
Iteration num (More iter, slower convergence, less stiffer)
The velocity update following projection is important to dynamic effects
Not only use in springs, but triangles, volumes, …
Pros & Cons:
(Usually in 3D softwares)
Pros:
Parallelable on GPUs (PhysX from NVIDIA)
Easy to implement
Fast in low res (Less than 1000 vertices)
Generic, can handle other coupling and constraints (including fluids)
Cons:
Not physically correct (no acc sol, stiffness related to meshes account and iter num)
Low performance in high res
Hierarchical approaches (From low to high res. But causes oscillation and other issues)
Acceleration approaches (Chebyshev, …)
Strain Limiting
Some normal updates + strain limiting (corrections)
e.g.: The constrain is not strictly equal to rest length but some constraints ( - stretching ratio as a limit)
Can be used to simulate some cloth whose stiffness increases when being stretched heavily.
The constraints can be used for reducing oscillation for FEM / …
Triangle Area Limit
Limit the triangle area. Define a scaling factor first (Mass center no change)
Properties:
Avoiding instability and artifacts due to large deformation
Useful for nonlinear effects (Biphasic way)
Help address the locking issue
Projective Dynamics
Main Idea
Instead of using projections to update directly as in PBD, projective dynamics uses projection to define a quadratic energy (identical to spring energy, but have 2 intermediate results & )
Finding the corresponding forces: (assumption: ) (identical to spring force)
=> althought the results are the same, the Hessians are different
The Hessian is the 2nd order derivative of energies to the vertices. The coefficients (of the diagonal) can be regarded as the numbers of the springs the vertices connected (3, 2, 3, 3 for vertex 0-3)
This matrix is simple and constant
=> Applying implicit integration
Simulation by PD
Combination fo the projectiive dynamics with implicit time integration. Use a direct LU solver and can factorize once (usually the most expensive part is the factorization and here the costs can reduce)
Principles
Newton’s Method calculate as Hessian
But in practice, usually not calculated exactly (The performance dep. on how well the Hessian is approx., identity matrix can even be used (but convergence takes longer))
Truncating the Hessian to make it p.d. (Lec. 5)
A (magical) diagonal approx (Lab 2)
A const matrix by PD (Lec. 6)
Pros and Cons
Pros
Have theoretical solution with physical meaning
Fast on CPUs with a direct solver (one factorization)
Fast convergence in the first few iter
Cons
Slow on GPUs
Slow convergence over time (fails to consider Hessian caused by projection) (suffering from high stiffness)
Cannot handle constraint changes easily (contacts / remeshing due to fracture / …)
Constrained Dynamics
This method is suitable for strong constraints / contacts / articulated rigid bodies (ragdoll animation) …
A critical problem: what if constraints/forces are very stiff (or even infinitely stiff)
example: multiple rigid bodies: human body can be regarded as several rigid bodies with very strong constraints (joints)
Compliant constraint:
The energy can be rewrited:
Let be the number of vertices and be the number of constraints ( - compliant matrix, - Jacobian, - Lagrangian multipliers (Dual variables))
By implicit Integration: (2 sets of var: primal var (or ) and the dual var )
Method 1: Solve by a direct solver together (primal-dual fashion): The matrix should be p.d. => need to make p.d.
Method 2: Reduce the system by Schur complement and solve first (not easy to reduce / not always sparse)
If reduce the first -> get the implicit method formula
Infinite stiffness in this case could be solved by (smaller makes the system better)
Lecture 7 Linear Finite Element Method (FEM)
Linear Finite Element Method (FEM)
The Linear FEM Assumption
In a nutshell, linear FEM assumes that for any point in the reference triangles, its deformed correspondence is (“Linear” for this linear formula)
For any vectors between 2 points we can use to convert it from ref to deformed state:
For example for edge 10 & 20: (indep on translation)
Problem: is related to deformation but contains rotations (want to also reduce rotation)
Green Strain
Ideally, we need a tensor to describe shape deformation only. Recall that SVD gives , where only and are relevant to deformation.
Get rid of as: (Green Strain, symmetric, only 3 var) ()
If no deformation: ; if deformation increases, increases
3 deformation modes:
is rotaion irrelavent (if adding rotaion then deformation grad is but is the same)
Strain Energy Density Function
Consider the energy density per ref area: :
Total energy:
The Saint Venant-Kirchhoff Model (StVK): ( and are Lamé parameters)
(Trace is the diagonal sum; - Second Piola-Kirchhoff stress tensor (some form of force density))
Force:
The are solved before, we concentrate more on the next part
Recall that:
By def
(the method is too complex => simplification)
For tetrahedron (3D ref -> 3D deformation)
Finite Volume Method (FVM)
For simple triangle / tetrahedra, FVM is equivalent to linear FEM
FVM considers force calculation in an integration perspective other than a differentiation perspective
Traction and Stress in FVM
Basic Ideas in FVM
General Case:
Traction: The internal force per unit length (area) => the total force at the interface is the integral:
Stress tensor : (projection: interface normal -> traction)
The Triangle Case: Take the midpoints which will be on the curve to estimate (the force contributes the same to the 3 vertices, the integral is for )
Force contributed by an element: (the integral of traction)
Divergence theorem: ( is constant within the elem)
Force: (Only is unknown)
In 3D (tetrahedra): FVM does the surfacial integral
Force:
Different Stresses
The stress tensor here: mapping from the interface normal to the traction (but kind of different from FEM)
In FEM: Stress , normal and traction are in the reference state
In FVM: this stress assumes the normal and the traction are in the deformed state
Output \ Input
Interface normal in the ref state (Undeformed)
Interface normal in the current state (Deformed)
Traction in the ref state (unformed)
Second Piola-Kirchhoff stress ()
-
Traction in the cur state (formed)
First Piola-Kirchhoff stress ()
Cauchy Stress ()
( - Deformation gradient)
The Finite Volume Method
Second Piola-Kirchhoff Stress: as in previous FEM formulation
=> To find (In the figure below, will be computed via )
The dot product of the edge 10 with
Meanwhile: & equal to 0
Algorithm (Explicit)
-> Various material (Heterogeneous materials) / various elastic / …
Hyperelastic Models
StVK: cannot handle flips / … (only widely used in graphics)
Hyperelastic: More used in engineering (material / mechanics /… ) => more general
Isotropic Materials
Isotropic
Force on any direction affects the same (with same behaviors) (-> ansotropic)
- Green Strain -> SVD; - Principal stretches: the singlar value of
In many lit, parameterize by principal invariants: ( is the right Cauchy-Green deformation tensor)
Models
The Saint Venant-Kirchhoff model (StVK):
The neo-Hookean model: (More standard model in material mechanics)
Two terms of : left - Against shearing; right - Against bulky change (volume change) (same as StVK)
The Mooney Rivlin model:
The Fung model:
Solving for First PK Stress ()
Use the principal stretch for computation => first P-K stress
Algorithm
Note: for the same model the first PK in this algorithm should equiv to
The Limitation of StVK
The compression (with high stress) and flip cannot be handled with StVK model (even get stable in the flipped side) => UNSTABLE
=> Neo Hookean solve the problem
Lecture 8 Linear Finite Element Method II
Nonlinear Optimization
Gradient Descent
Solving => The negative gradient dir is the fastest descent dir
Algorithm
Step Size Adjustment
Linear search methods: Exact line search (solve another optimization problem) & Backtracking line search ( is a smaller num (~ 0.3-0.4) -> Wolfe Condition)
Descent Directions
The dir is descending if a sufficiently small step size exists for: (at least descending)
ie. (the same side with the negative dir)
With line search. use any search dir as long as it’s descending
Descent Methods
Gradient descent:
Newton’s method (if the Hessian is always positive definite):
Any method using a positive definite matrix to modify the gradient yields a descent method:
A unified descent framework:
Differences between these methods:
Total cost = Per-iter cost * Num of iter
Lecture 9 Collision Handling
Collision Detection
Collision Detection Pipeline
Broad-Phase: Remove the pairs which are impossible to collide
SH: Easy to implement; GPU friendly; needs to recompute after updating obj
BVH: More involved; not GPU friendly; update BV to update BVH
Narrow-Phase: Which pairs really collide
Board-Phase Collision Culling
Spatial Hashing (Partitioning)
Static space division
Spatial partitioning divides the space by a grid and stores objects into grid cells
Cell 12:
Cell 13:
Cell 14:
Cell 15:
Cell 8:
Cell 9:
Cell 10: ,
Cell 11: N/A
Cell 4: ,
Cell 5: ,
Cell 6:
Cell 7:
Cell 0:
Cell 1: N/A
Cell 2:
Cell 3:
For example for : Cell 5, 6, 9, 10 => other triangles stored in these cells: & => pairs: & and do further detection
For moving objects, just expand the object region (add more)
Optimization for memory costs: Memory wastes in empty cells and too many meshes should be built in 3D. Instead of allocating memories to cell, build an object-cell list and sort them
One question is how to define the cell ID. Using the grid order is not optimal, since it cannot be easily extended and it is lack of locality (LHS). Morton code uses a Z-pattern instead (RHS).
Bounding Volume Hierarchy (BVH)
Bounding volume hierarchy is built on geometric/topological proximity of objects.
External Object
The parts that are not like to collide in different boxes => test if box intersection occurs
To find elem potentially in collision with an obj, just traverse the tree:
Internal Intersections (Self Collision)
2 procedures in the tree
xxxxxxxxxx
Process_Node(A)
{// check if the parental node (A) has self-collision
Two AABBs intersect if and only if they intersect in every axis; shperes: need to compute the distance
Narrow-Phase Collision Test
Discrete Collision Detection (DCD)
Actually not collision detection but intersection detection
DCD tests if any intersection exists in each state at discrete time instant:
To a triangle mesh, the basic test is edge-triangle intersection test: easy and robust
(The means the interpolation value between the points other than time)
Tunneling
The difference between collision and intersection detection:
Intersection detection only considers if the 2 triangles (/…) have intersection at the two discrete time steps
Collision detection can consider the penentration
Problem: The green triangle at state and doesn’t have intersection with the blue one but actually the position relationship changes => tunneling problem (penetrating through each other)
Not usually occurs for bulky objects but can happen for thin surfaces
Continuous Collision Detection (CCD)
CCD tests if any intersection exists between two states: and
To a triangle mesh, there two basic tests: vertex-triangle and edge-edge tests
Vertex-Triangle:
( which is a function about time)
Need to solve a third order equation (with one variable) => has formulae solution (not recommended, need to solve cubic root) / use binary division (0, 1) and get the first sol
Edge-Edge:
Issues with CCD:
Floating-point errors, especially due to root finding of a cubic equation
Buffering epsilons, but that causes false positives
Gaming GPUs often use single floating-point precision
Computational costs: more expensive than DCD
Some argue that broad-phase collision culling is the bottleneck
Difficulty in implementation
Continuous Collision Response Approaches
The Two Continuous Collision Response Approaches
Given the calculated next state , we want to update it into , the path from to is intersection-free
Interior Point Method: Use the rest state to the collided state to estimate the every next state (safer)
Slow: far from solution; all vertices; cautiously by small step size
Always succeed
Impact Zone Method: Use the final point to estimate, then do optimization in the “impact zone”
Fast: Close to solution; only vertices in collision (impact zones); can take large step sizes
May not succeed
Interior Point Methods (Cont.)
Log-Barrier Interior Point Methods
For simplicity, just consider the Log-barrier repulsion between 2 vertices (usually in graphics do a cutoff at some distance)
Implementation
Formulate the problem:
Gradient Descent:
The step size must be adjusted to ensure that no collision happens on the way. To find , we need collision tests.
Impact Zone Optimization (Cont.)
The goal of impact zone optimization is to optimize until it becomes intersection-free. (This potentially suffers from the tunneling issue, uncommon)
Such that:
The optimizations are always trying to get close to the solution
Rigid Impact Zones
The rigid impact zone method simply freezes vertices in collision from moving in their pre-collision state. It’s simple and safe, but has noticeable artifacts
Only used when all other methods not work
A Practical System
Untangling Cloth (Discrete)
Intersection Elimination
Consider how to eliminate existing intersections, but without using any collision history
Useful when there are already intersections in simulation, due to:
Past collision handling failures
Intense user interaction
In this case, we don’t require the simulation is to always intersection-free
Eliminating cloth-volume and volume-volume intersections is straightforward: simply pushing vertices/edges in the volume out
Untangling Cloth
Cloth-cloth intersection is complicated (don’t have a clear definition of inside and outside)
Baraff et al. used flood-fill to segment cloth into regions and decided which region is in intersection. (Cannot handle boundary well, not friendly in GPU)
divide into several regions => which is larger => the smaller part moves to the larger side
Volino and Magnenat-Thalmann proposed to untangle cloth by reducing the intersection contour
Their method can handle boundaries, but it doesn’t always work
Lecture 10 Surface Waves
Two Types of Simulation Approaches
Lagrangian Approach: Dynamic particles or mesh (Node movement carries physical quantities (mass, velocity, …))
In 2D, a (1.5D) height field is a height function (cannot have 2 at one position) and the velocity is also a function (negative -> opposite direction)
The height field and velocity can be represented by: ( means the volume flowing through a position at the unit time)
Ignoring advection and external acceleration, we get a simple form of velocity derivative: ( - density; - pressure)
From to : => the velocity changes due to the change of pressure (also dep on density)
increases when and increases when
Shallow Wave Equation
Simplification of the two equations by assuming shallow wave (the wave is small and the height change can be neglected)
Eliminate and formulate the shallow wave equation (indep on the velocity field)
Finite Differencing
Discretization
Finite Difference
Use the difference to approximate the derivative
Forward differencing (first order, 2nd order error)
Backward differencing (first order)
Central differencing (second order)
Second-Order Derivatives (central differencing in the fluid field)
Applying central differencing twice to estimate height (Laplacian)
Also for pressure
Discretized Shallow Wave Equation
Want to obtain the height field in the next time step
Problem: Volume loss or generation
Volume Preservation
We want the volume to stay the same
Suppose , but by using the formuation in the prev section, we will get
and we cannot guarantee is 0
Solution 1
Modify scheme into: (not using as the coefficient => using and instead)
Can be understanded as the water exchanges between and . The water flowing from the LHS will flow into the RHS => Volume preserved
Solution 2
Assuming in the right term is const
Pressure
The pressure is related to the water height:
Replace the term with a constant
Viscosity
Viscosity tends to slow down the waves (Add a viscosity constant : too high: not viscos; too low: not like fluid)
Boundary Conditions
Dirichlet Boundary: The boundary height is constant (open boundary)
Neumann Boundary: specifies the boundary derivatives. For example, a zero-derivative boundary (no water exchange) means (closed boundary)
Algorithm with Neumann Boundary
– to 3D –
Two-way Coupling
The coupling will be two ways: between liquid->solid and solid->liquid (also can be bubbles / soft solids / …)
Key: how to expel water out of the gray cell regions
Virtual Height
Idea: Set up a virtual height so that
Poisson’s Equation
The outcome is Poisson’s equation, with and being unknowns
Regard as a Laplacian operator for every equation ( & are all 0, no virtual height)
Algorithm
Process: Set boundary condition -> Mask (which grids need solving) -> PCG_Solve() to retrun -> Update the height field (Usually need to multiply a factor, prevent too large waves when dragging the cube <- explicit’s unstability)
Rigid Body Update
Estimate the floating (buoyancy) force by the actual water expelled in every column
Also need to consider rotation => torque
Lecture 11 Eulerian Fluids
A Grid Representation and Finite Differencing
A Regular Grid Representation
Storing physical quantities over a grid => Fields
Scalars: Density/color, Pressure, Temperature, …
Vectors: Velocities
Finite Differencing on Grid
Grid is good for central differencing
(Same for the direction)
Discretized Laplacian
Boundary Conditions
The boundary condition specifies if it’s outside
Dirichlet:
Neumann:
Example: Laplace’s Equation
With the grid, discretize Laplace’s equation: (every grid has a function , every )
Note: at least one condition should be Dirichlet condition (avoid inifinite solutions)
Laplacian Smoothing is basically average the self value with neighbors => the distribution will be smoother and smoother
(for simulation, time can be added along the term)
The process of applying Laplacian smoothing is called diffusion
Staggered Grid
Problem with Central Differencing
Central differencing gives the derivative in the middle
The cell doesn’t exist at
To get we need and => weird because is unused
Solution: Staggered Grid
We define some physical quantities on faces, specifically velocities
The x-part of the velocity is defined on vertical faces
The y-part of the velocity is defined on horizonal faces
They represent the flow speed between two cells, at cell :
Divergence-Free Condition
No volume change is equal to say the fluid is incompressible (divergence-free velocity field)
Bilinear Interpolation
Use bilinear interpolation to interpolate physical quantities (3D: trilinear)
Method of Characteristics: solving a long PDE in steps
Step 1: Update by solving
Step 2: Update by solving
Step 3: Update by solving
Step 4: Update by solving
Step 1: External Acceleration
Straightforward, just add acceleration to and
For gravity as ext. acc.: ( for verticle velocities <- affected by gravity force)
Step 2: Advection
Solving in an Eulerian way can be unstable => use the real meaning of advection: carry physical quantities by velocity (this problem does not exist in Lagrangian methods) => Semi-Lagrangian Method
Semi-Lagrangian Method
Trace a virtual particle backward over time
Steps: (same for (verticle))
Define
Compute
Compute
Note: if the velocities are staggered, we need to do staggered bilinear interpolation
=> Subdivide the time step for better tracing
Steps: (Make => do 4 times of updates in an original timestep)
Define
Compute
Compute
Compute
Compute
Step 3: Diffusion
The laplacian of velocity ( and )
Note: if is large: can be unstable => use smaller sub-steps / implicit integration
Step 4: Pressure Projection
The pressure is caused by incompressibility => By updating the pressure, we should achieve divergence free condition:
Eventually got a Poisson equation:
With bc:
Dirichlet (open):
Neumann (close):
Once solve => update and done
Air and Smoke
Air Simulation
Steps: (can be used in underwater as well)
Step 1: Update flow (velocity field)
Step 2: Use semi-Lagrangian advect all of the other physical quantites (density / temperature / …)
Typically use Dirichlet boundaries for an open space and Neumann for a container
Water Simulation
Water simulation also involves air-water interface. Water doesn’t occupy the whole space => Need to solve the air-water boundary
Representations:
Volume-of-fluid (VOF) (every grid has a percentage of fluid and air) => not accurate, not tells the interface
A signed distance function defined over grid
Advect: (needs corrections)
Semi-Lagrangian (volume loss)
Level set method (spec. to signed distance function update) (volume loss)
Volume loss affects significantly when perturbation is large